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We take further steps in the development of the characteristic approach to enable handling the 
physical problem of a compact self-gravitating object, such as a neutron star, in close orbit around 
a black hole. We examine different options for setting the initial data for this problem and, in order 
to shed light on their physical relevance, we carry out short time evolution of this data. To this end 
we express the matter part of the characteristic gravity code so that the hydrodynamics are in con- 
servation form. The resulting gravity plus matter relativity code provides a starting point for more 
refined future efforts at longer term evolution. In the present work we find that, independently of the 
details of the initial gravitational data, the system quickly fiushes out spurious gravitational radia- 
tion and relaxes to a quasi-equilibrium state with an approximate helical symmetry corresponding 
to the circular orbit of the star. 

PACS numbers: 04.25.Dm, 04.20.Ex, 04.30.Db, 95.30.Lz 

I. INTRODUCTION 

The problem of computing the evolution of a self-gravitating object, such as a neutron star, in close orbit about 
a black hole is of clear astrophysical importance, both to understand systems thought to drive such spectacular 
phenomena as gamma-ray bursts (see, for instance, 0) and to predict the details of the gravitational radiation which 
could be observed by the new generation of gravitational wave detectors (see, for instance, Furthermore, the 

dynamics of the system could be quite rich. For instance, the tidal interaction between the black hole and the 
star could be strong enough to tidally disrupt the star with a consequent drastic change in the emitted waves 
Alternatively, as has recently been suggested Q, only a portion of the star's mass might be transferred to the black 
hole and angular momentum transfer might boost what remains of the star into a wider orbit. The system would 
then undergo another inspiral phase and this process might repeat itself. Consequently, the expected gravitational 
waves would display a very rich structure of several 'chirp-phases'. 

A better understanding of these possibilities requires, as basic building blocks, solving Einstein equations coupled to 
an appropriate matter field without restrictive assumptions. To this end, long-term well behaved numerical simulations 
must be available. Presently, several numerical relativity codes for treating the two body problem are either being 
developed or planned [ESS 111 j as the know-how for simulating Einstein equations becomes more mature. Although 
most of these efforts concentrate on the Cauchy approach to Einstein equations, there are alternatives approaches 
which have been successful for specific problems. In particular, the characteristic formulation of general relativity 
has shown remarkable robustness to deal with single black hole space-times. Stable axisymmetric studies of Einstein 
equations coupled to perfect fluids have been achieved 0, 0; and our previous work has produced stable three- 
dimensional characteristic numerical relativity codes for vacuum space-times 'll^ and for fluid space-times with very 
small pressure 12]. 

This paper develops further methodology necessary for the evolution of a self-gravitating star or other object in 
orbit near a Schwarzschild black hole. Towards this final goal, we examine the issue of setting initial data. In either 
the characteristic or Cauchy approaches to this problem, a serious source of physical ambiguity is the presence of 
spurious gravitational radiation in the initial gravitational data. Because the characteristic approach is based upon a 
retarded time foliation, the resulting spurious outgoing waves can be computed by carrying out a short time evolution. 
We carry out such a study in the present work. We find that, independently of the details of the initial gravitational 
data, such spurious waves quickly radiate away, and that the system relaxes to a quasi-equilibrium state with an 
approximate helical symmetry corresponding to the circular orbit of the star. This result provides important physical 
justification of recent app roaches for initializing the Cauchy problem which are based on imposing an initial helical 
symmetry [3 El El lH El . 

We also examine two useful tools which can be applied to long term simulations as well as to the initial data 
problem: (i) a tool to monitor the development of a singularity in the null coordinates and (ii) a tool to use co- 
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rotating coordinates so that an orbiting star remains approximately at a fixed coordinate position. 

A crucial ingredient of any hydrodynamical simulation is the proper handling of shocks and discontinuities and 
the proper conservation of baryonic mass. Robust numerical techniques are available to this end, which express the 
system in a conservation form designed for handling discontinuities in the fluid Ts'l . In particular, high resolution shock 
capturing schemes utilize the fluid's characteristic propagation speeds to capture the discontinuities in an accurate 
way. In the present work, we follow the formalism presented in |19| |. where the general relativistic hydrodynamic 
equations are presented in a conservation form adapted to both the Cauchy and the characteristic formulations. High 
resolution shock capturing schemes have been successfully incorporated in the characteristic approach in the case 
of axisymmetric space times 0, 0| . Although this is the ideal approach to the hydrodynamic problem, here we 
implement a less time consuming algorithm due to Davis po| . This approach has limitations in dealing with shock 
formation but is sufficient for our present purpose to study the initial data problem. 

The initial data problem for the characteristic formulation of relativity Tl'l has received little attention in comparison 
with that for the Cauchy problem. This is partly because the data for the characteristic formulation is unconstrained, 
whereas the data for the Cauchy formulation has to satisfy an elliptic system of constraints. Even so, in both 
formulations, care must be taken to set matter data that represents the intended physical problem. In particular, it 
is not trivial to set gravitational data free of spurious gravitational radiation. A method for initializing characteristic 
data based upon a correspondence with Newtonian theory 0,11^ guarantees that the resulting gravitational radiation 
obeys the Einstein quadrupole formula in the Newtonian regime |23l l23 |. However, as discussed in the domain 
of applicability of the method is limited to when (i) the Newtonian potential $ of the star satisfies 4> <C 1 and (ii) 
relative velocities are much slower than the speed of light. The second issue could be easily dealt with by making the 
orbital radius a of the star large compared to the black hole mass M. However, the following practical reasons make 
it desirable to initialize the star in a close orbit with a < lOM (for which v > c/3): 

• Since the orbital period of the star around the black hole scales as a^^^, a simulation starting from large a 
would require very time consuming runs in a regime that could be described approximately by less expensive 
perturbative treatments. 

• When using spherical coordinates, the resolution at large a worsens when using a uniform grid. Computer 
resources limit the number of angular grid points that can be used. Thus, adequate resolution of the star 
requires that a be small. 

• The most interesting physics, intractable by other means, occurs when the compact object is at small a. 

Construction of mathematically consistent initial data is much simpler than construction of physically meaningful 
data. While the former involves the solution of constraint equations (if any), the latter is less mathematically explicit. 
Not only does one want matter data describing an orbiting star but also gravitational data with minimal spurious 
radiation. Success must be gauged through actual evolution of the data and analysis of the resulting space-time. [4^ 

There are basically no constraint equations for the gravitational initial data in the characteristic formulation so 
that mathematically consistent data is trivial to prescribe. In order to investigate the physical problem we proceed 
by investigating gravitational data based upon two completely different underlying assumptions: 

• Gravitational initial data obtained by means of the Newtonian correspondence method 

• Spherically symmetric gravitational data corresponding to a shear-free initial null hypersurface. 

The first method applies only to a quasi-Newtonian regime. The second method ignores the focusing effect of the 
star which introduces shear in its nearby null rays. Nevertheless, in both cases the gravitational field relaxes to a 
state of approximate helical symmetry after about a light-crossing time for the system. These results support the 
view that the details of the initial gravitational data may not be important as long as they are reasonable and some 
time is evolved in which the spurious radiation flushes out. 

The plan of the paper is as follows: in Sec. ^ with present a summary of previous work on the characteristic 
formulation of numerical relativity and on conservative hydrodynamics. The construction of initial matter and grav- 
itational data is covered in Sec IIIII in both stationary and co-rotating frames. In Sec. IIVI we discuss the conditions 
under which coordinate singularities might develop. Some details of the numerical implementation are presented in 
Sec. Computational tests and results are given in Sec. IVII 
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II. SUMMARY OF PREVIOUS RESULTS AND NOTATION 



A. Characteristic formulation of Einstein equations 

The formalism for the numerical evolution of Einstein's equations, in null cone coordinates, is well known [Tlll2fill2^ 
(see also [26, 27, 28]), and is based upon the analytic treatments of 29, 30, 31]. For the sake of completeness, we give 
here a summary of the formalism, including some of the necessary equations. The version of the gravity code being 
used here is fully described in [s^, 113 ■ It has most recently been applied in jsj . 

We use coordinates based upon a family of outgoing null hypersurfaces. We let u label these hypersurfaces, 
{A = 2,3), label the null rays and r be a surface area coordinate. In the resulting = {u,r,x^) coordinates, the 
metric takes the Bondi-Sachs form I291I30I 



ds^ = -{e^^il + —)-r''hABU^U'']du^-2e^^dudr-2r^hABU''dudx^ 



r 

+ r'^hABdx^dx^ , (1) 

where h^^hsc = ^"^^ det{hAB) ~ det{qAB), with qab a unit sphere metric. We work in stereographic coordinates 
x^ = {q,p) for which the unit sphere metric is 

QABdx'^dx^ = 4^{dq^ + dp^), where F:^l + q^+p^. (2) 

Our previous work used P — 1 + + p^, here we change notation to F because we will use P to represent pressure. 
We also introduce a complex dyad q"^ — -|-(l,i) with i — For an arbitrary Bondi-Sachs metric, Hab can then 

be represented by its dyad component 

J = hABq\''/'2, (3) 

with the spherically symmetric case characterized by J = 0. We introduce the (complex differential) eth operators 3 
and § 113, as well as a number of auxihary variables K = hABQ'^q^ U = U^qA, Qa = r^e~^ /^HabU^ , Q — Qa^^, 
B = ^(3,v = m and k = dK. 

The Einstein equations Gat = decompose into hypersurface equations, evolution equations and conservation laws. 
Here we note that the hypersurface equations form a hierarchical set for z/,., fc^^, /3,r, B ^, Q.r, U.r and W^r, and the 
evolution equation is an expression for {rJ)^ur- The explicit form of the equations is given in |32l ] in the vacuum case 
and in |0 for the matter terms. Note the following correction to the matter source term in Eq. p^-f31') 



2 {rJ).^^,^ - ((1 + r-^W) {rJ)^^ = -r'^ [r'^W) ^ + 2r-^e'^?i^e^ - [r-^W) ^J + Nj 

{{JVang " KVangf + KL) , (4) 



where Vang — vaQ"^ |12| (with Va the velocity) and Nj are defined in [Tlj and |32| . The remaining Einstein equations 
are the conservation conditions which are satisfied here because of the simple choice of boundary data. 

The null cone problem is normally formulated in the region of space-time between a timelike or null worldtube F 
and 2^. We represent 2^ on a finite grid by using a compactified radial coordinate x — r/{l + r). The numerical 
grid is regular in {x,q,p), and consists of two stereographic patches covering the north and south hemispheres, each 
containing rix x Uq x np grid points. The a:— grid covers the range [0.5, 1], and each angular grid patch extends at 
least two grid points beyond the equator so that there is an overlap region. 

The mass of the Schwarzschild black hole is denoted as M, and in all the computational tests we will take M — 1. 
The star has mass m and radius i?*. 



B. Characteristic hydrodynamics in conservation form 



The integration of the hydrodynamical equations is done in a more accurate way if the system can be expressed 
in conservation form. This allows for better conservation of baryonic mass and the possibility of exploiting the 
characteristic structure of the equations to resolve shocks via Godunov methods p^ . Here we use the formalism 
developed in jl^], where the equations Va J° = 0, WaT'^^ = are expressed in conservation form as, 

a.o iV^v^) + a,. iV^F^^) = , (5) 
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where = {D,S'^,E) are "conservative" variables defined as, 



and the fluxes and sources are 



D 


= UO 


= J" = 
















E 


= 


_ 2^00 


= - 





(6) 
(7) 
(8) 



pw^ 



= J' 

= T^' = phu'u^ + Pg'^ , 
pji = ^ phu°u^ + Pg°\ (9) 

S" = 0, 

S'' = -V^Tl^T^\ (10) 

The state of the system is uniquely described in terms of the geometry, the primitive matter variables {p, P, h, u°) or the 
conservative variables, the fluid's equation of state and the normalization condition u^Ua = — 1- In the particular case 
of a perfect fluid which we consider here, the relation between primitive and conservative variables is straightforward 
and does not require expensive inversion methods (see |l9l | for details). 



III. INITIAL AND BOUNDARY DATA 



There exist only rough physical guidelines for prescribing initial matter and gravitational field data in general 
relativity for a star orbiting a black hole. Here we present initial data which is at least mathematically consistent with 
Einstein's equations and has some underlying connection with the Newtonian picture of a star in equilibrium, which 
is orbiting a collapsing central object. The Newtonian picture is not a good approximation to the relativistic regime 
in which we run our simulations. Nevertheless, the justification for this approach is that the evolutions, presented in 
Sec. lVIl relax on the order of a light crossing time to a more astrophysically realistic state, as extraneous gravitational 
wave content in the initial data is radiated away, and the gravitational field adapts to the approximate symmetry of 
the matter distribution. This relaxed state is remarkably independent of the initial gravitational data. 



A. The initial matter data 



We prescribe the initial matter data within a simple Newtonian framework. We take the star to be a spherically 
symmetric polytrope of index n = 1 |2SJ with 



for i? < i?*, and p = for i? > i?*, where p is the density at radius R. Denoting the pressure by P, the equation of 
state is 



m sin .§2L 

R. (11) 



P=^. (12) 

TT 

Note that the maximum value of P/ p is Tn/2R^, so that Newtonian theory gives a good approximation to an equilibrium 
configuration provided the polytrope is not near its Schwarzschild radius. We prescribe the initial matter velocity 
to be uniform across the polytrope. In the evolutions considered later the center of the polytrope will be placed at 
{q = p = 0,r = a), with the velocity set for a circular orbit, so that = = initially. 

The assumption that the polytrope is spherically symmetric is always valid mathematically since the initial data 
is freely specifiable, but the star remains in equilibrium only if there is no tidal force. Of course, it is known how to 
compute the tidal distortion of a poly trope in equilibrium within both Newtonian gravity and post-Newtonian general 
relativity - see for example [stL l38j. However, in the simulations of a star orbiting a black hole presented in this 
paper, the polytrope is not stable against tidal disruption, i.e. the simple Roche indicator e/j = R^j a^J^M jm 3> 1. 
Thus, for the present work, there is no point in calculating an equilibrium configuration for the tidal distortion. 
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In order to obtain physical characteristic initial data for the problem we transform the above Newtonian initial data 
into density, velocity and gravitational fields within the framework of a Bondi-Sachs metric in general relativity. As 
discussed in Sec. this needs to be done in a relativistic regime in which the Newtonian correspondence method 0, 
123 . is only roughly approximate. 

The density and velocity matter fields are regarded as being given in a Lorentzian frame in the Minkowski space- 
time in which the center of mass of the polytrope is at rest. Then a coordinate transformation is made from this 
local Lorentz frame to the global Bondi-Sachs coordinates. This approach is strictly valid only if space-time curvature 
can be neglected in the vicinity of the star, which requires that the radius of the star must be small and that its 
gravitational source effect must be negligible. The approach has a proper Newtonian correspondence in the limiting 
case of small velocities, small m/i?» and large a. 



1. The density and velocity fields 

We prescribe initial data (p, V"') for the localized distribution of matter described in Sec. IIII Al The matter is 
described in a (locally) inertial frame Sj = {t,x,y, z) in which the center of mass of the matter is instantaneously at 
rest at the origin. Globally, we have Bondi-Sachs coordinates centered about the black hole (with mass A/). To the 
extent that the gravitational effect of the matter can be neglected, the geometry is spherically symmetric and can be 
described in Eddington-Finkelstein coordinates Se = (u,'', (JiP)- The origin of Si is a,t u — 0,r = a,q = p = and, 
for simplicity, we choose a central matter velocity of the form (T^",0, y*,0) in Se- Except in the case M = 0, the 
curvatures of Si and Se are different, and thus it is not possible to construct a unique global transformation between 
Si and Se- However, we can construct a transformation that is valid near the origin of Si. 

We proceed by constructing a locally inertial frame S' = {t' , x' , y' , z') whose origin is at O = {u — 0,r ~ a,q ~ 
0,p = 0) and with the S' axes pointing in the {u,r,q,p) directions. In general. Si is moving relative to S'. When the 
black hole mass M = 0, the space-time is flat and there is an unambiguous transformation from 5" to Se 

t' 
I 

X 



li + (r — 


a) 


2rq 




1 + g2 4- 




2rp 




l + q^^ 




2r 




1 + g2 4- 


- p^ 



(13) 



When M 7^ 0, the metric is Schwarzschild in the Eddington-Finkelstein form 



,2 



ds' = -[l \du' ~2dudr + - —{dq'+dp')- (14) 

V r / (1 + 9 +P ) 

An orthonormal tetrad aligned with S' and defined in Se at the point O is given by 



Xl,^ = ((l-f^)"^0,0,0) 



Xl) = (0,0,0,^-) 



2M\-5 r 2M 
a J 

In a neighborhood of O, Se and S' are related by the transformation 



- (-(l--) ,0,0). (15) 



2M\-\ , / 2M\ 

t'- 1 

a / V a / 

2M\\ , 

^ 

a / 



X 

y' 

^ = 2^' 



(16) 
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from which we construct the inverse transformation 

2Af\5 / 2AfN 



t' = 1 "+ 1 kr-a) 



a J V a / 

x' = 2aq 
y' = '^ap 

z' = (l-— ) {r-a). (17) 

In the M — case, the transformation is unambiguous and gfobal, but when M ^ Q the transformation to a focaUy 
inertial system can only be defined in a neighborhood of O. Nevertheless we need a transformation that is valid in a 
region around O, which we construct by combining and H17|) to obtain 

2M\k / 2M\-i 



t' = 1 u + 1 {r-a) 



r J 



I 

X 



2rq 



1 + g2 + p2 
2rp 



2r \f 2M\-h 
^— ^-r-a 1 . (18) 

The transformation H18|l reduces to H13|l when M = 0, and to H17|) near the point O. 

It is most convenient to describe the matter in an inertial frame in which it is at rest which, in general, is not the 
case for 5". In order to achieve this we need to make a Lorentz transformation to coordinates Si — {t, x, y, z). In the 
case that the matter velocity is entirely in the g-direction, which includes the case of a quasi-circular orbit around the 
black hole, the transformation between Se and Si is 

t = j{t' - vx'), X = -/{x' - vt'), y^y\ z = z', (19) 

where v is the velocity between Si and S", and 7 is the usual Lorentz factor, 7 = (1 — v'^)~^ . Then combining 119|) 
and lfT7|l we find 

2M\-i, , 2rqv \ 

r / 1 + + / 

"('-T)"'<"->) 




(20) 

The transformation H20|l is sufficient for setting scalar initial matter data, such as the density. Given an Se grid at 
M = 0, we find for each point the coordinates {x,y,z). The density is then determined by its values in Si. We also 
need to set the initial 4-velocity. In Si the 4-velocity is = (—1,0,0,0). Thus in Se the covariant 4-velocity is 

Vim = ^yil) ^_ dt 



Q^(E)a b Q^(E)a 
2M\k 



2qv / 2Mn^-5 ^ (r - a)M 2M\^ 



(■ 



_ 3 



1 + g^+p^ V r ) V r J 

2rv{l+p'^ - q^) 

(1 + g2 +p2)2 

Avqpr 



(1 + g2 +p: 



,2-12 



(21) 
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In order to set a value of v for a quasi-circular orbit, we need V^i^"* and l/f^)"^ at O. We find 

2M\i / 2M\-i 



A circular orbit has zero radial acceleration. Thus, applying the geodesic equation, we find 



M a-2M 
„r^'^-V^^ (24) 

As expected, a circular orbit can exist only if a > 3M. 

B. Initial gravitational data 

There is not a well-developed theory for initializing the gravitational field of a star in close orbit about a black hole. 
We consider two options here. The first is to simply set 

J(u = 0, r, x^) = 0. (25) 

In the absence of matter, this choice would eliminate all spurious gravitational waves but in the presence of matter it 
is physically unrealistic. It implies that the null rays generating the initial null hypersurface are shear-free, i.e. there 
is no bending of light by the star as would expected in a normal astrophysical scenario. This introduces spurious 
gravitational waves which superimpose with the gravitational field of the star to cancel the bending effect. The second 
choice of J is quasi-Newtonian gravitational data 21, 22, 23, 2^] determined by a Newtonian correspondence method 
which introduces the correct bending effect and radiation content in the Newtonian limit. 

This quasi-Newtonian data is astrophysically realistic only when the Newtonian potential and the matter velocity 
are small. These conditions are only very roughly satisfied for the relativistic binary considered here. The Newtonian 
potential approaches unity as r approaches 2M and the velocity of the body approaches that of light. Nevertheless, it 
is reasonably simple to compute quasi-Newtonian initial data and it is interesting to compare the resulting evolutions 
with those initialized by J = 0, e.g. to compare their spurious radiation content and the quasi-equilibrium states to 
which they relax. 

The procedure for determining quasi-Newtonian data involves solving a sequence of Poisson equations which, in the 
Newtonian limit, lead to exact agreement with the Einstein quadrupole formula for the initial gravitational radiation 
content in the system. Here we apply the method only at leading order where Eq. (3.8) of |23l| gives (in the present 
notation) 



where $ is the Newtonian potential satisfying 



{r^J^r).r ^ -29^$ (26) 



= Attp (27) 



in terms of the Euclidean Laplacian obtained in the Newtonian limit. The computation of $ is straightforward and 
Eq. (|26(l can then be solved by a combination of numerical and analytic techniques subject to boundary conditions 
that determine the integration constants. Equation H26(l establishes a Newtonian correspondence on the initial null 
hypersurface u = 0. The sequence of successive Poisson equations ensure that (I26() is satisfied in time in terms of a 
Taylor expansion in u. 

At the r = 2M surface of the black hole, we choose the boundary conditions 

J J,, = at r = 2M (28) 



for solving the Poisson equations. For Af = 0, these boundary conditions reduce to regularity conditions at the vertex 
of the outgoing null cones. Details of the calculation are given in the Appendix. 
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C. Inner boundary data 



Our aim is to simulate a star which is in orbit close to an astrophysical object undergoing gravitational collapse. We 
model the exterior field of the collapsing object by the vacuum geometry of the outgoing u = const null hypersurfaces 
whose inner boundary is a white hole horizon. In terms of the Kruskal picture of an isolated Schwarzschild black hole 
this inner boundary would be the past branch of the r = 2M hypersurface. Thus our simulations correspond to the 
characteristic evolution of data on two intersecting null hypersurfaces, the initial outgoing null hypersurface Aq and 
the white hole horizon 7i~. Since the orbit of the star is exterior to 7i~ the inner boundary data for the matter is 
trivial. 

The construction of general gravitational data for such a double-null initial value problem has been reduced to the 
integration of propagation equations (ODE's) along the null generators of 7i~ [2^. The free data on 7i~ consist of 
the specification of its conformal metric, represented by J, and its shift represented by the Bondi variable U . 

Here we make the simple choice that J = vanishes on 7i ~ . One of the propagation equations (the Raychaudhuri 
equation) then implies that the intrinsic expansion of 'H~ is constant along its generators. We choose the initial value 
of this expansion to vanish so that r = const on 7i~. We then set r — 2M on the sphere S~ where and Aq 
intersect. Thus the inner boundary 7i~ is described by 



These choices simplify the integration of the remaining propagation equations on 7i~. Choosing the "shift" of 
to vanish leads to the further simplification that C/ = on 7i~. We consider this choice first. (Below we consider the 
shift corresponding to coordinates co-rotating with the orbiting star.) The remaining data at S~ is the initial twist 
of the horizon, which determines U^r on the horizon and the initial outward expansion of S~ , which determines the 
initial value of /3. In the spirit of simplicity, we set the twist to zero, which assigns zero angular momentum to the 
white hole, and we set (3 to zero, which in the pure Schwarzschild case sets the scale of retarded time u on the horizon 
to the standard inertial time measured by observers at null infinity. 

Because all free initial inner boundary data at Sq have been chosen to be identical to the inner boundary data 
of a Schwarzschild horizon in standard Eddington-Finkelstein form, and because we have chosen J = on 7i^, this 
data propagates along the generators of to give the Eddington-Finkelstein boundary values for the corresponding 
Bondi variables, namely J = U = U^r = f3 = Q,W — —2M. Alternatively, we could prescribe Minkowskian inner 
boundary data by setting the value of the inward expansion to match that of a ingoing (collapsing) Minkowski light 
cone. 

The Bondi evolution system also requires the value of J.^ on 7i~. However, for the simple boundary data presented 
above, the horizon propagation equations imply that — const on 7i~, and our boundary condition (|28|l then 
implies J^. = 0. 



It is also desirable to be able to specify data a.tr — 2M such that the angular coordinates of a star in uniform circular 
orbit around the black hole with angular velocity Q remain constant. In terms of a standard angular coordinate 0, 
this can be arranged by the transformation 9—^9 + Qu on the data presented above. This introduces a shift on 
the horizon which is represented b y a non-zero value of C/. At large distances from the horizon this shift becomes 
superluminal, but calibration tests [40j have shown that this does not preclude solving the characteristic initial value 
problem. In the Schwarzschild case, the vector field du in the non-rotating coordinates is the time translation Killing 
vector T, whereas in the rotating coordinates du equals the helical Killing vector T + fl^ where $ is a rotational 
Killing vector. 

The calculation of the required U, which must be carried out in the stereographic (q^p) angular coordinates used in 
the code, was performed by means a computer algebra script. We start with the Bondi-Sachs metric in standard form 
for a Minkowskian space-time, with coordinates {r,q,p,u), and transform {r,q,p) to Cartesian coordinates {x,y,z), 
as specified in Eq. H13() (here the Cartesian coordinates are written unprimed, even though they are written primed 
in Eq. We then perform a rigid rotation about the y-axis through an angle 9, leading to Cartesian coordinates 

(a;', y', z') with 



J at r = 2M. 



(29) 



1. 



Co-rotating data 



x' = X cos 9 + z sin 9 



y ^ y 



z 



— 2: cos 
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— X sm 



61, 



(30) 
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Finally, we transform back to Bondi-Sachs coordinates {r' , q' , p' , u') by computing the transformation {r,q,p,u) 
(r' , q' ,p' ,u') and the Jacobian for the transformation when is small; then these are employed to find the metric 
of Minkowskian space-time in {r' ,q' ,p' ,u') coordinates. In order to check the calculation, via a computer algebra 
script we checked that all components of the Riemann tensor of the transformed metric are zero. Next, we apply 
the coordinate transformation {r,q,p,u) — > ir' ,q' ,p' ,u') to the Se metric H14|) . We find that the metric remains in 
Bondi-Sachs form with J = (3 — Q,W — —2 and 

2fr\l + q^ ~p^) 
" (l + g2+p2)2 

(1 -I- q^ ^-p-^Y 

where / = d9/du, and for convenience the ' symbols have been dropped. Thus, 



^_^1±^!^V±|^. (32) 
1 + q-^ +p^ 



The value of / is set by the condition that the center of the polytrope should be at rest with respect to the ' coordinates. 
We achieve this by transforming the covariant velocity H21|) to the ' coordinates, and then evaluating the contravariant 
velocity at r — a, q = p = 0. Wc find 



V^^ ^=[fa-vJl^'^-^]; (33) 



thus = if we set 



When V is given by Eq. (|24|l . 



then Eq. (|32|l becomes 



1 _ 2M 



(34) 



f^\H- (35) 



U^-M'-±f^P^. (36) 

V a3 1 + gf2 + p2 ^ ' 

Note that U^r = 0. Furthermore, the transformation to ' coordinates does not change the formulas for Vr, Vq and Vp 
given in Eq. 1)21(1 : the formula for Vu is changed, but that is not part of the required initial data. 



IV. REGULARITY OF THE NULL COORDINATES 



The Bondi description of an asymptotically flat space is based upon the surface area coordinate r, which is assumed 
to increase monotonically to infinity along the outgoing null rays. This monotonicity is built into the initial data by 
assuming that the matter variables and the gravitational data J are well-behaved functions of r, for 2M < r < oo. 
However, certain choices of matter data can lead to gravitational data that is unrealistic astrophysically. For example, 
consider a spherical star of mass m and radius i?* located at r = a, with a » 2M so that the linear approximation 
for the bending of light by the star is valid. Then, under normal astrophysical conditions, when 



the outgoing null rays which graze the surface of the star would be bent by the star so that they would travel to infinity 
along asymptotically parallel trajectories. This lack of expansion of the outgoing rays would be a breakdown of the 
regularity of the r-coordinate, and it would also introduce a large shear in the geometry of the outgoing null rays. 
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This loss of regularity can be mathematically avoided by prescribing gravitational data with zero (or small) shear 
along with the matter data for the star. But clearly that would not represent an astrophysically realistic problem. 
Note that, for a given star mass m, this breakdown occurs for large values of a and not when the star is close to 
r = 2M. 

It is important to monitor this effect. An afhne parameter A measured along the radially outgoing null rays must of 
course increase monotonically in an asymptotically flat space. Thus the relevant quantity to measure is the expansion 
d\r^ which in a Bondi coordinate system is given by 

d),r = (38) 
In our setup, d\r = 1, i.e /3 = 0, at the inner r = 2M boundary, but the Bondi hypersurface equation 

Or 13 = ^{J,rJ,r " K^) + 2TTr{p + P)v^ (39) 

8 ' 

implies that (3 increases outward monotonically. Thus the expansion dxr decreases monotonically and reaches its 
smallest value at infinity. Initial data for which i9a'"|oo ~ 1/10, corresponding to Poo ~ 1, would either represent an 
extreme astrophysical scenario (such as a star about to enter a black hole) or would signal an imminent breakdown 
of the r-coordinate. 

In Sec. IVI El we illustrate how this effect can be monitored by measuring f3oo- This is important in order to 
avoid wasting computational time trying to simulate systems which are either unrealistic astrophysically, or where a 
coordinate singularity will develop. 

V. NUMERICAL IMPLEMENTATION 
A. Hydrodynamical equations 

We implement an algorithm of Davis 20], which can be regarded as adding artificial dissipation where needed to 
a MacCormack scheme, by means of a slope delimiter procedure. To do this requires knowledge of the largest of the 
eigenspeeds (speed of characteristic modes of the principal part of the fluid's evolution equations) in each direction. 
The expressions for these can be read-off from those worked out in . 

This dissipation, shown in ID for simplicity, takes the form 

where U"'^^ is the update from the MacCormack step and -Di+1/2 is defined as. 



A+i 



with 



/2 = (^t+i/2(^^) + ^7+1/2(^7+1)) (t^JVi - u7) (41) 



(42) 



j+1/2 



<Ac/;ii/2'^f^;ii/2> 

ic(z;)(l-*(r±^,/j) . (43) 



The values oi v, C and 5* are 



V = max{ I A j I } — 



dt 
dx 



Civ) 



*(r) = 



v{l~v) if t; < 0.5 
0.25 otherwise 



min(2r, 1) if r > 
otherwise 
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B. Update strategy 

In this work we combine the PITT characteristic vacuum code developed in [Til I25I [s^ (and thoroughly tested and 
applied in a variety of situations, see for instance 0,|3E3j0|) with the general relativistic hydrodynamic equations 
provided above. Although this is the first application of the combined equations in three-dimensional settings, we 
follow closely the strategy pursued successfully in two-dimensional scenarios 0, ^3 • The hierarchy of integration of 
the equations is basically the following: 

1. With data on an initial hypersurface Afu, the metric is updated to the new level J^u+Au- Here the fluid variables 
at the intermediate level A/'„_|_^„/2 1 needed for the integration of the metric equations, are approximated, to first 
order in time, by their values at A/'«. Note that since the typical propagation speeds of the fluid is less than the 
speed of light, this approximation in general is acceptable. 

2. Next, the general relativistic hydrodynamic equations are updated to Afu+Au- 

Although the integrations of the GR/hydrodynamical equations are written in second order form when the fluid/GR 
variables are frozen, the above procedure is formally only a first order accurate approximation. Formally higher order 
schemes can be obtained by iterating several times per time-step, where "intermediate" values of the fields are 
employed (consisting of the average of field values at Mu and the previously obtained fields at Afu+Au- Alternatively, 
a more efficient scheme can be devised by keeping an extra level of the fluid variables (at Mu-Au): so that one can 
extrapolate to second order their values at Afu+Au/2- In the present work, since we concentrate in rather short time 
evolutions, we have opted not to do this in order to reduce the amount of time required by the code. Additionally, 
during the time evolutions we use a F-law equation of state given hy p — (T — l)pe, with F = 1 -I- 1/n and ignore the 
effects of viscosity and magnetic fields, since their dynamical time scales are much longer than those considered here. 

VI. COMPUTATIONAL TESTS 

Our main goal here is to study the initial evolution phase of a star in orbit around a black hole as a preliminary step 
toward carrying out long term evolutions. Of primary concern is whether the system can be initialized in a way which 
allows a meaningful simulation of the ensuing inspiral and capture of the companion star. The strategy is to begin 
with a rough choice of gravitational data, and then show that the system quickly relaxes to a state which provides 
physically reasonable matter and gravitational data, which can in turn be used to initialize a longer evolution. 

Several options are available for carrying out this study in a discriminating and efficient way. They involve the 
choice of initial data; the choice of coordinates fixed with respect to the black hole or co-rotating with the orbiting 
star; the possibility of setting the black hole mass to zero; and the use of switches in the code that allow running the 
gravitational field or the hydrodynamics in a frozen mode (which expedites the turn-around of the numerical tests). 
In all tests, the star is initialized with a uniform angular velocity about the black hole corresponding to the circular 
orbital velocity of its center. We then choose from the following specific options: 

• Initial gravitational data given by either J ~ or quasi-Newtonian data. 

• Full evolution or evolution with the internal hydrodynamics of the star frozen, i.e. the gravitational field reacts 
to a rigidly orbiting source. 

• Fixed coordinates (C/ = at the inner boundary) or co-rotating coordinates (U given by the value calculated in 
Sec. lmTTTll . 

• Inner boundary consisting of either a mass M = 1 Schwarzschild event horizon or of an ingoing null cone in 
Minkowski space-time (no black hole). 

We also carry out experiments to check the range of stellar mass, size and orbital radius for which light bending 
effects do not introduce singularities in the null coordinates. The observed behavior of the code is described in the 
sub-sections below. Unless otherwise indicated, all norms refer to the L2 norm. 

A. Code convergence 

The convergence of the PITT code has been established in a series of papers for evolution in vacuum or with fluids 
having negligible pressure 0, IT^ Is^ . l40l | . The new ingredient here is the addition of a fluid obeying a polytropic 
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Q 




FIG. 1: For the case m = 10"^ = 0, i?* = 3, a = 9 and grid sizes given by 45^ x 63, 65^ x 93 and 85^ x 123, the diflFerent 
self-convergence factors are obtained. The behavior observed is consistent with at least first order convergence. 

equation of state. The fluid is treated numerically by a first order accurate algorithm. To illustrate this we place the 
star in a space-time without a black hole and monitor the solution under different grid resolutions. 

Because the evolution of the hydrodynamics with the present serial code is considerably time consuming, the 
convergence test must be based upon limited grid resolution. We evolved the case m = 10~^, i?* = 3, a = 9 using 
quasi-Newtonian initial gravitational data. Three different grid sizes were used: 45^ x 63; 65^ x 93 and 85^ x 123 
(stereographic x radial). With this choice of grid sizes, the middle and higher resolution grid spacings correspond 
to 2/3 and 1/2 the grid spacing of the base grid, respectively. In order to assess the convergence behavior of the 
implementation we examine the following quantities. 



Qf3 



Q global 



|/3(A) - |/3(A2/3)| 
|/3(A2/3)-|/3(A/2)| 

. f |f,(A)-|f,(A2/3)| \ 
"^^n|F,(A2/3)-|f^z(A/2)|/ 

/ S,,-(fz,-(A)-fz,-(A2/3)f 
Vs/,(i^/,(A2/3)-F,,(A/2)f 



(44) 
(45) 

(46) 



The capital index / ranges over all gravitational fields, i.e. Fj = J, U, W}, and the index j over all points in the 
base grid which are common to the middle and high resolution grid. The factors Q^, Qmin and Q global measure the 
convergence of the field /3, the minimum of the convergence rate obtained for all gravitational field quantities Qmin 
and a global value Qgiobai including all fields, respectively. The results are displayed in figure^ where the horizontal 
lines provide the value for first and second order convergence. The observed behavior indicates a convergence rate 
consistent with first order convergence. 



B. Evolution for different gravitational initial data 



Independently of the initial data J, the star settles into "quasi-equilibrium" after a couple of hydrodynamical times 
given by 2i?*, i.e. all variables relax to roughly constant values with respect to a co-rotating observer. Of the small 
remaining variation, a major portion subsequently dies off after a couple of crossing times 2a, the time needed for the 
star to communicate with the inner boundary and back. This latter behavior is observed whether or not the internal 
hydrodynamics of the star is evolved. After this relaxation period, the time dependence of the gravitational field is 
mainly determined by the motion of the star. 

As an illustration of this behavior we evolve a star in a space-time both with and without a black hole. The star 
is placed at a = 9 with mass m = 10~^ and radius R^, = 3 and the gravitational initial data is given either by 
J — (shear-free data) or by the quasi-Newtonian value. For purposes of comparison, runs are made both with 
the hydrodynamic variables evolved and frozen. Figures |21 and |31 illustrate the behavior of J for these cases. After 
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FIG. 2: Behavior of || J|| for the case with no black hole. Three scenarios arc plotted: Full evolution (hydrodynamics + gravity) 
using initial data J — (dashed line with "x" symbols); full evolution using quasi-Newtonian initial data (solid line with 
circular symbols); and gravitational evolution with frozen internal hydrodynamics using quasi- Newtonian data (dotted lines 
with triangular symbols). After some transient behavior, lasting until about w ~ 2, all curves roughly approach the same value. 
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FIG. 3: Behavior of || J|| for the case of an M = 1 black hole. Again three scenarios arc plotted: Full evolution using initial data 
J = (dashed line with "x" symbols); full evolution using quasi-Newtonian initial data (solid line with circular symbols); and 
gravitational evolution with frozen hydrodynamics using quasi-Newtonian data (dotted lines with triangular symbols). After 
some transient behavior lasting until w ~ 2, all curves roughly approach the same value. 



a relatively brief transient behavior, the norm of J approaches a value quite independent of the initial data (with 
similar behavior observed in all other gravitational field variables). This indicates that most of the spurious radiation 
present in the initial hypersurface is 'flushed out' in a short time. In order to elucidate whether the system settles 
into a quasi-stationary state, as might be expected, we analyze next whether there is an approximate helical Killing 
vector in the space-time. 
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FIG. 4: L2 norm of the difFerence between the computed values of J for the quasi-Newtonian and spherically symmetric initial 
data. 

Agreement of the "relaxed state" 

As indicated above, the system appears to relax to a state independent of the details of the initial data. In order to 
quantify this, we calculate the L2 norm of the difference between the two numerical solutions obtained with Newtonian 
and with shear-free initial data. Figures ^ and [S] illustrate the behavior of the difference vs. time, showing how as 
time progresses the agreement between the numerical solutions becomes more pronounced. 

It is also illustrative to observe the point- wise difference between the two numerical solutions. In particular, 
inspection of \JNewt — Jspher\ at X"*" reveals not only that the solutions tend to agree as time progresses, but also that 
this agreement is more marked along null rays passing through the star and neighboring rays. 

C. Quasi-equilibrium behavior 

In order to measure the approach to quasi-equilibrium we monitor the rate of change of the angular part of the 
metric along the streamlines of the vector field ^ = T -I- kH.^. As explained in Sec. IIII (TTl for k = 1, ^ equals the 
helical Killing vector of the background black hole with set to the initial orbital angular velocity of the star. For 
comparison purposes we also consider k = 0, for which ^ equals the static Killing vector T of the background black 
hole, and k = —I, for which ^ equals the helical Killing vector counter-rotating with respect to the orbital motion of 
the star. We measure the change of the gravitational field with respect to the fiow of these vector fields by the norm 
Fi^ = WC^hABW^ ■ For orbital motion in quasi-equilibrium around the black hole, we should then find -F'k=i ~ 0. In 
order to verify that this is the case, we monitor the values of for k = 1, 0, —1. Again, the star is initialized at a = 9 
with mass m = 10^^ and radius i?* = 3. 

The first test, carried out with quasi-Newtonian initial data and shown in Figure El compares the behavior of 
for the three different values of k indicated above. At first, the curves differ very little, but at later times there 
is a marked difference (> 1.5 orders of magnitude), indicating that the system approaches an approximate helical 
symmetry. For the k = I helical case, Fi decays more than two orders of magnitude from its initial value. 

Next we repeat the test with J = shear-free initial data, with the results shown in Fig. [7| The initial values of 
are not the same as those in the first test, but a similar time behavior is observed. Again the helical choice Fi decays 
more than two orders of magnitude from its initial value, and more than an order of magnitude more than either Fq 
or F-i. 
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FIG. 5: Absolute value of the difference between the computed values of J at X"*" evolved from quasi-Newtonian and from 
shear-free initial data. The panels (A,B,C,D) correspond to the time sequence u — 0, u = 0.75M, u = 1.5M, and u — 2.5M. 
Panel A has been truncated for comparison purposes, as its height is fa 9 x 10~^. 
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FIG. 6; Comparison of F^, using quasi-Newtonian initial data. 

D. Co-rotating coordinates 

Co-rotating coordinates can improve the tracking and hydrodynamic treatment of the orbiting star by reducing 
or eliminating the variation of field values due to purely coordinate effects. They are introduced here by specifying 
the value of U at the inner boundary, as described in section IlII C II By defining U via Eq. H36|l . we have checked 
that this indeed keeps the angular coordinate of the star fixed. The tests were carried out in the star mass range 
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FIG. 7: Comparison of Fk. using shear-free initial data. 

m G [10"'^, 10"^], with M = 1, i?^ = 3M and a = 9M. Figure |H1 illustrates the results for m = 10~^. As can be 
seen clearly from the density contours displayed, the co-rotating coordinates maintain the central density of the star 
at the same initial coordinate location after evolution through u = 2.8. Note that, although the initial polytrope 
is spherically symmetric about the center of the star, the kinetic energy contributed by the initially uniform orbital 
angular velocity skews the density distribution. These results verify that the initialization can be performed equally 
well in co-rotating coordinates. 

E. Monitoring the expansion of the null coordinates 

When the star is initialized sufficiently far from the central black hole, its bending effect on neighboring light rays 
can cause a coordinate singularity, by reducing the expansion of the outgoing null hypersurfaces to zero. As discussed 
in Sec. IIVI the bending effect on outgoing null rays is first evident at future null infinity, where it is manifested in 
our formalism by /3 ^ oo. However, /? — > oo also when the generators of the outgoing null hypersurfaces approach a 
black hole horizon, in which case /3 ^ oo on a complete spherical set of rays. As a result, it is difficult to distinguish 
at a given retarded time between whether it is a horizon or a coordinate singularity that is responsible for the lack of 
expansion. In the present context, what matters is whether such coordinate singularities can arise on a relatively short 
dynamical time scale which would have bearing on our main conclusion that the system relaxes to a quasi-stationary 
orbit. To assess this possibility, we carry out evolutions with a star of mass m = IQ-^M and radius = 3M 
initially placed in orbit at different distances a from the central black hole. We then monitor the expansion e^'^^ for 
roughly the "relaxation" time. Figure |51 graphs the w-dependence of the minimum value of e~'^^ over the sphere at 
null infinity for the initial orbital radii a = 9M, 16. 5M, 20M, 25M, 30M. As expected, the expansion becomes smaller 
as the separation between the star and the black hole increases. More important to the initialization problem, for 
a given separation the minimum expansion does not appear to change considerably from its initial value during the 
relaxation time. An estimate obtained with the lens equation combined with the initial value of f3 at infinity seems 
to provide a good indicator of whether coordinate singularities will develop. 

We now examine the possibility of dealing with a more massive star m = lO^^M, keeping i?» = 3M. The lens 
estimate predicts zero expansion at a separation given by a > 22. 5M. To study this problem, we set up data at 
different separations and, in view of the behavior just seen in the case with m = 10~^M, we evolve for a short 
time (u = M/10). Figure [TUI graphs the computed value of e~'^^ vs. a, for a range of separations. As the star is 
placed farther from the black hole, the expansion gets considerably smaller. Furthermore, the values obtained for 
a > 16AI only represent an upper bound since even with the largest practical grid nr x nq x np = 132 x 82^ the 
star is considerably unresolved. As a result, the neighboring null rays pass by the star at a significant "grid distance" 
and do not accurately represent the expansion that would be calculated in the continuum problem. Nevertheless, the 
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FIG. 8: For the case ni — 10^ ', i?, = 3M, a — 9M we show density profiles in the plane with stereographic coordinate p — Q. 
The motion takes place in the {r,q) plane. In all figures, the vertical axis is x (the compactified radial coordinate) and the 
horizontal axis is q. The top panels correspond to the fixed coordinates (A) and co-rotating coordinates (B) at the initial time 
M = 0. Figures (C) and (D) show the corresponding density profiles at w = 2.8. The co-rotating coordinates, panel (D), perform 
well in keeping the coordinate location of the star fixed. For comparison, panel (C) shows the actual displacement of the star. 

results in the figure suggest that it would be possible to simulate a star in orbit around a black hole with a < 16M 
on a uniform grid without developing coordinate problems. 

F. Code performance 

We illustrate the speed of the code for the case of the finest grid used in the convergence test, i.e. 81^ x 123. Setting 
the time-step to Au = 0.014M, a run until u = 1.5M takes 9 hours on a 2.4GHz Pentium 4 processor, and requires 
1.4Gb of memory. Prom this we conclude that a single orbit would take roughly 11/2 months. 



VII. CONCLUSION 



Within the characteristic framework, we have developed and implemented a numerical relativity code, as well as 
procedures for finding the required initial data, for evolving a star in close orbit around a Schwarzschild black hole. 
Wc have shown that after a short evolution time the system relaxes to a quasi-equilibrium state which is mainly 
independent of the initial gravitational data. Variations of the initial matter data, such as the shape or size of the 
star, have not been investigated. For small variations about stellar equilibrium one would expect similar results since 
the accompanying variations in the star's gravitational field should relax to the same quasi-equilibrium state. We 
have also developed and demonstrated tools that permit the use of co-rotating coordinates and that monitor possible 
problems with the null coordinate system. 
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FIG. 9; Time behavior of the minimum expansion e " at future null infinity vs. initial location of the star. As the star is 
placed farther from the black hole its focusing effect on the central null rays increases. However, the expansion stays close to 
its initial value for a given separation a. 
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FIG. 10: Behavior of the minimum of the expansion e at future null infinity vs. location of the star. As the star's initial 
location is placed farther from the black hole its focusing effect on the central null rays becomes considerably. 



Successful simulation of a neutron star in close orbit around a black hole requires a gravity code, a hydrodynamic 
code and physically ap prop riate initial data. It has already been demonstrated that the characteristic gravity code 
is accurate and stable |lll l32| | . In this paper we have implemented and tested a conservative formulation of the 
hydrodynamics and used it to shed light on the problem of prescribing physically meaningful initial data. There are 
several issues that remain to be addressed before physically interesting long term evolutions can be carried out. The 
results obtained in this paper would appear to justify the effort required to make these improvements. 

Specifically, in order for the code to yield astrophysically useful results, three further conditions must be fulfilled. 
First, more realistic inner boundary data must be provided. This must correspond to a spinning black hole and 
must include the gravitational distortion induced by the companion star. Second, the turn-around of results must 
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be considerably improved. This will require revising the numerical algorithms and parallelization of the code to take 
advantage of large platforms. Finally, since accurate simulations require that the numerical error be well below the 
expected radiation output (< 5%), the use of adaptive mesh refinement seems crucial. Preliminary work in this 
direction has recently been presented of these items are major tasks and are deferred to future work. 
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APPENDIX A: COMPUTATION OF THE QUASI-NEWTONIAN INITIAL DATA 



The Newtonian potential outside the polytrope is simply 



Inside the polytrope, the derivation is rather more complicated. We find 



Rtt 



{R < i?,)- 



and from it obtained 3^$. First, we make the following definitions 

z — q + ip 



Pr, 



{r — of' + zz{r + aY 
1 + zz 



Then we may write 



<3 



= -12t3 {R>R,) 



a^r^z^m^/P^ 



t 



5/2 



1 — zz 



{--K^r^ + 2r7r^a h 3i?^ - 7r^a^)sint2 - 3i?*7r, / — cost; 



P. 



(R < R^ 



(Al) 



(A2) 



(A3) 



(A4) 



Finally, J is obtained by numerically integrating (r'^J^r),r — — 2S'^$ in second order form starting from the inner 
boundary r = 2M. 
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